The accuracy of integrated rotational quaternion and angular velocity from curvature

The problem of integrating rotations and angular velocity from a known timeand space-dependent curvature is met in the dynamic analysis of non-linear spatial beams, if the numerical formulation is based on curvature, as in [1]. To avoid singularities and make the formulation more convenient, rotational quaternions q̂ are used for the parametrization of rotations. The relation between the curvature quaternion, κ̂ , and the rotational quaternion in the quaternion body basis is shown in Eq. (1): κ̂ = 2q̂ ◦ q̂. (1) Here q̂∗ is the conjugated rotational quaternion, and the prime ( ′ ) denotes the space derivative. For a given curvature, Eq. (1) constitutes the homogeneous system of four differential equations with variable coefficients and, together with the initial condition for q̂, forms an initial-value problem for q̂. This is essentially the same problem as integrating rotations from the angular velocity, where ω̂ = 2q̂ ◦ ̇̂ q. (2) Here ω̂ is the angular velocity quaternion, and the dot over the symbol denotes the derivative with respect to time. The problem of integrating rotations from the angular velocity is met in, e.g. navigation, robotics and computer graphics, and several numerical solution methods have been proposed [2, 3, 4, 5]. However, an exact closed-form solution is only possible for the angular velocity being constant in time (or the curvature constant in space). The majority of present numerical methods fail to preserve exactly the unity constraint of the rotational quaternion, and experience a sudden instability due to the horizontal stretching of quaternion response curves. The identical rotational quaternion must be obtained from integrating either the given curvature or the related given angular velocity, which implies that the two quantities are not independent. Their mutual relation reads: ω̂ ′ = 1 2 (ω̂ ◦ κ̂ − κ̂ ◦ ω̂)+ ̇̂ κ. (3) Here ω̂ ′ is the space derivative of the angular velocity and ̇̂ κ is the time derivative of the curvature. Scalar parts of the angular velocity and curvature quaternions are equal to zero, which is why Eq. (3) represents a system of three nonhomogeneous differential equations with variable coefficients. When the curvature is a given function, Eq. (3) along with the initial condition for ω̂ constitutes an initial-value problem for ω̂. If the angular velocity is given, Eq. (3) represents an initial-value problem for κ̂ . Only if certain conditions are satisfied, exact solutions for such systems of differential equations with variable coefficients can be found [6]. For convenience Eqs. (1)–(3) are written in the component form and rearranged into matrices having varying, yet known coefficients, and arrays of unknown functions. The solution for the rotational quaternion is presented in Eq. (4), where q is an array of four unknown components of the rotational quaternion, and K4D is a given 4x4 skew-symmetric matrix: q(x2, t) = e 1 2 ∫ x2 x1 K4D(ξ q(x1, t), q0 = q(x1, t). (4) The solution for the angular velocity is presented in Eq. (5): ω(x2, t) = e − ∫ x2 x1 K3D(ξ ,t)dξ ω(x1, t)+ ∫ x2 x1 e ∫ x2 η K3D(ξ ,t)dξ k̇3D(η , t)dη , ω0 = ω(x1, t). (5) Here ω is an array of three unknown components of the angular velocity, K3D is a known 3x3 skewsymmetric matrix, and k̇3D is the time derivative of its axial vector. The condition for the solutions to be exact [6] becomes rather simple if we consider that the matrices K4D and K3D are skew-symmetric. It only requires that the components of the curvature be arbitrarily scaled functions of the same form. It can also be proven that even if this condition is not met, the solution approaches the exact solution when the observed interval [x1,x2] decreases. Hence, with the use of an appropriately small interval, Eqs. (4)–(5) represent approximate solutions, which converge toward the exact solutions in the limit. In order to test the accuracy of the present approximate solutions, two rotational quaternions are chosen in an analytical form as functions of both x and t, from which the exact curvature and angular velocity are obtained using Eqs. (1) and (2). The present numerical solutions are compared to the analytical ones on the domain {(x, t): 0 6 x 6 10; 06 t 6 5 }. Small computational steps ∆x = x2− x1 and ∆t = t2− t1 are used. In the first example, the quaternion is chosen such that the curvature satisfies the condition for the solution to be exact. In this case, the present and the analytical solutions match within the machine round-off error precision for all components of both the rotational and the angular velocity quaternions, regardless of the size of the computational step. The unity constraint of the rotational quaternion is perfectly preserved. The curvature in the second example does not satisfy the exactness condition. The unit norm of the rotational quaternion is nevertheless preserved for any length of the integration step. The results for both the rotational quaternion and the angular velocity converge monotonically if ∆x is decreased, with no horizontal stretching taking place, as shown in Fig. 1. Decreasing ∆t does not make errors lower, because all integrations in Eqs. (4)–(5) are conducted with respect to x. The maximum normalized errors of the rotational quaternion and the angular velocity are of the same order of magnitude.